{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np  # 导入包\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50  #生成50个样本点\n",
    "x = np.linspace(0, 10, nsample)  # 从0-10之间生成50个数\n",
    "X = np.column_stack((x, x**2))  # 生成第一列是1，第二列是x，第三列是x的次方\n",
    "X = sm.add_constant(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.72063062, 2.098422  , 2.98932883])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta = np.array([5, 2, 3])\n",
    "e = np.random.normal(size=nsample)\n",
    "y = np.dot(X, beta) + e\n",
    "model = sm.OLS(y, X)\n",
    "results = model.fit()\n",
    "results.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   1.000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   1.000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>1.875e+05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 10 Nov 2020</td> <th>  Prob (F-statistic):</th> <td>2.00e-92</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>21:38:35</td>     <th>  Log-Likelihood:    </th> <td> -75.079</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    50</td>      <th>  AIC:               </th> <td>   156.2</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    47</td>      <th>  BIC:               </th> <td>   161.9</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    4.7206</td> <td>    0.457</td> <td>   10.332</td> <td> 0.000</td> <td>    3.801</td> <td>    5.640</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    2.0984</td> <td>    0.211</td> <td>    9.931</td> <td> 0.000</td> <td>    1.673</td> <td>    2.524</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>    <td>    2.9893</td> <td>    0.020</td> <td>  146.287</td> <td> 0.000</td> <td>    2.948</td> <td>    3.030</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.967</td> <th>  Durbin-Watson:     </th> <td>   1.919</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.374</td> <th>  Jarque-Bera (JB):  </th> <td>   1.894</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.409</td> <th>  Prob(JB):          </th> <td>   0.388</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.510</td> <th>  Cond. No.          </th> <td>    142.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       1.000\n",
       "Model:                            OLS   Adj. R-squared:                  1.000\n",
       "Method:                 Least Squares   F-statistic:                 1.875e+05\n",
       "Date:                Tue, 10 Nov 2020   Prob (F-statistic):           2.00e-92\n",
       "Time:                        21:38:35   Log-Likelihood:                -75.079\n",
       "No. Observations:                  50   AIC:                             156.2\n",
       "Df Residuals:                      47   BIC:                             161.9\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          4.7206      0.457     10.332      0.000       3.801       5.640\n",
       "x1             2.0984      0.211      9.931      0.000       1.673       2.524\n",
       "x2             2.9893      0.020    146.287      0.000       2.948       3.030\n",
       "==============================================================================\n",
       "Omnibus:                        1.967   Durbin-Watson:                   1.919\n",
       "Prob(Omnibus):                  0.374   Jarque-Bera (JB):                1.894\n",
       "Skew:                           0.409   Prob(JB):                        0.388\n",
       "Kurtosis:                       2.510   Cond. No.                         142.\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**R方值已经非常精确了**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAFlCAYAAAA+t0u5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de1jUZf7/8efNMCqlRh4qwWNmhkWeqDC3NG2zzJSs1LLTbr91D53bLM36dtrCos2yLcu143qujOywUWqWaZboqKhEqRkymFlJaqIwM/fvDwYWFQM5fWaG1+O6uGBuPjPzZq5tX9735z4Yay0iIiISGqKcLkBERET+R8EsIiISQhTMIiIiIUTBLCIiEkIUzCIiIiFEwSwiIhJCop0uAKBVq1a2Y8eOTpchIiJSb1auXPmjtbb1we0hEcwdO3YkMzPT6TJERETqjTHmu4raNZQtIiISQhTMIiIiIUTBLCIiEkJC4h5zRYqLi8nLy2Pfvn1Ol1KnmjRpQtu2bXG73U6XIiIiISBkgzkvL49mzZrRsWNHjDFOl1MnrLX89NNP5OXl0alTJ6fLERGREBCyQ9n79u2jZcuWERvKAMYYWrZsGfGjAiIiUnUhG8xARIdyqYbwN4qISNWFdDCHmgceeIAnnnjisL9PT09nw4YN9ViRiIhEmogJ5nSPl74TF9Fp3Hv0nbiIdI+3/mtQMIuISA1FRDCne7yMn5eFt6AQC3gLChk/L6tWwvmRRx6ha9eunH/++eTk5ADw73//mzPOOIPu3btz2WWXsXfvXpYtW8b8+fMZO3YsPXr0YNOmTRVeJyIi8lsiIpjTMnIoLPYf0FZY7CctI6dGr7ty5Upmz56Nx+Nh3rx5rFixAoDhw4ezYsUK1qxZQ0JCAi+++CJnn302Q4cOJS0tjdWrV9O5c+cKrxMRkTDz+eeQmlryvR6E7HKpI5FfUHhE7VW1ZMkSLr30Uo466igAhg4dCsC6deu49957KSgoYM+ePQwaNKjC51f1OhERCU2fvDKfs8aMINpXjC/azRdT59Lv+qF1+p4R0WOOi405ovYjUdGs6euvv55//etfZGVlcf/99x92uVNVrxMRkdCT7vGy8rW3aFS8n2gbINpXzMrX3qrzOUwREcxjB3Ulxu06oC3G7WLsoK41et1zzz2Xt956i8LCQnbv3s0777wDwO7du2nTpg3FxcXMmDGj7PpmzZqxe/fusseHu05EREJfWkYOpmg/UUAAQ7Erms/iT63xbdLKRMRQdkrPeKDkQ8wvKCQuNoaxg7qWtVdXr169GDlyJD169KBDhw6cc845ADz88MOcddZZdOjQgcTExLIwHjVqFH/605+YPHkyb7zxxmGvExGR0JdfUEiXH7eyoVVH3ks4h887nM6q+ARMDW+TVsZYa+v0DaoiKSnJHnwec3Z2NgkJCQ5VVL8a0t8qIhIu+k5cRP7OX2m2fy+7mjQta4+PjWHpuAE1fn1jzEprbdLB7RExlC0iIlKrpk/nvh7NaNLIfUAo18Zt0soomEVERMpbuhSuu44L336J1OGJxMfGYCjpKacOT6zxbdLKRMQ9ZhERkVrxyy8wejR06ACPPUZK8+Z1HsQHUzCLiIiU+tvfIC8PliyB5s0dKUHBLCIiDVa6x1u2oufq3C94eNZMePBB6NPHsZoUzCIi0iCVnrNQuqVz+vGn0arftXS86DqGOViXJn9VIi8vj2HDhtGlSxc6d+7MrbfeSlFREYsXL2bIkCGHXP/uu+/Ss2dPunfvTrdu3XjhhRccqFpERCpTes5C0tb13Lx0Fl1+zGVS8ggeX7jJ0boUzL/BWsvw4cNJSUnhm2++4euvv2bPnj1MmDChwuuLi4sZM2YM77zzDmvWrMHj8dC/f//6LVpERKokv6CQXt5sZs2+hzs+m8GM2RPo5c2u8TkLNRVZwVzLJ4AsWrSIJk2a8Ic//AEAl8vFpEmTeOmllyo8wnH37t34fD5atmwJQOPGjenatW7Xu4mISPXExcZwxdoFRAf8GMDt95Gcm1Ur5yzURPjcY66o5zliRMkMur17oW9fWLsWAgGIioLTT4dbb4Xrr4cff4TLLz/wuYsXV/qW69evp3fv3ge0NW/enPbt27Nx48ZDrm/RogVDhw6lQ4cODBw4kCFDhnDllVcSFRVZ//4REYkE95x1HGc9vBwL+E0Uxa5oPCf2qPMNRCoTPsFcmV9+KQllKPn+yy81fklrbYWnSx2uHWDatGlkZWWxYMECnnjiCT766CNeeeWVGtciIiK1yFoufvpeAkV7Sb1sLO4d29nYLYmRY4bX+7rlg4VPMP9WD/eoo2DGDBg4EIqKoFGjksel091btapSD/lgp556Km+++eYBbbt27WLr1q107tz5sM9LTEwkMTGRa665hk6dOimYRURCzc6dsGULUY9NZMIddzhdzQEiZ4y1Tx9YuBAefrjkey2sQRs4cCB79+7ltddeA8Dv9/P3v/+d66+/nqOOOuqQ6/fs2cPicv8AWL16NR06dKhxHSIiUstatIAvvoDbbnO6kkNETjBDSRiPH19rC8ONMbz11lu8/vrrdOnShZNPPpkmTZrw6KOPArBw4ULatm1b9uXxeHj88cfp2rUrPXr04P7771dvWUQklPz6K9x1V8ntzsaNS+YkhZjwGcp2SLt27XjnnXcOae/fvz+FhYdOqS89s1lERELQbbfBiy/C4MEVTyoOAaH3TwUREZG6MHcuTJsG48aFbChDFYLZGNPEGPOlMWaNMWa9MebBYHsnY8wXxphvjDFzjDGNgu2Ng483Bn/fsW7/BBERkUps2QJjxkBycsle2CGsKj3m/cAAa213oAdwoTEmGXgMmGSt7QLsBG4IXn8DsNNaexIwKXidiIhIvUv3eBnz12fY2Pt37N27nw8nPAlut9Nl/aZKg9mW2BN86A5+WWAA8Eaw/VUgJfjzsOBjgr8faA636Lfy967O08JKQ/gbRUSckO7xMmfyXJ6edicdd+bjCvh5+c3lpHu8Tpf2m6p0j9kY4zLGrAZ+AD4CNgEF1lpf8JI8oHRFdjywFSD4+1+AlhW85hhjTKYxJnPHjh2HvGeTJk346aefIjq4rLX89NNPNGnSxOlSREQiTlpGDv03fIbb7yPaWlwBPz03ryYtI8fp0n5TlWZlW2v9QA9jTCzwFpBQ0WXB7xX1jg9JV2vtVGAqQFJS0iG/b9u2LXl5eVQU2pGkSZMmtG3b1ukyREQiTlFePiPWfgiAL7jl5vL2iY4fUlGZI1ouZa0tMMYsBpKBWGNMdLBX3BbID16WB7QD8owx0cAxwM9HWpjb7aZTp05H+jQRERHw+5nywZM08RVzx8W303bXDpa3T2RVfALxDh9SUZlKg9kY0xooDoZyDHA+JRO6PgYuB2YD1wFvB58yP/j48+DvF9lIHo8WEZHQ88gjJG3yMGHIbcw/9byy5hi3y/FDKipTlXvMbYCPjTFrgRXAR9bad4G7gTuMMRspuYf8YvD6F4GWwfY7gHG1X7aIiMhhLF5csiRq9GjOePDvxMfGYID42BhShyc6fkhFZSrtMVtr1wI9K2jfDJxZQfs+4IpaqU5ERORINWoE550HU6aQ0qwZKb3Cax6PtuQUEZHIcvbZsGCB01VUm7bkFBGRyPD44yUHVPj9TldSIwpmEREJf0uXwj33lGy9GYInRh0JDWWLiEhYSvd4ScvIYe+27fz31dto1qYtR//731C9zSZDRnj/s0JERBqkdI+X8fOyOH7dKuZOv4uWu37iugvuIH3znsqfHOIUzCIiEnbSMnJI2LKOmbPv4aSf88BAoKgo5LfbrAoFs4iIhJ38gkKSc7OIDvgxgLGW5NyskN9usyoUzCIiEna6Re/HWEuxK/qAfbDjQny7zarQ5C8REQkvfj8vLZxM7Ipl3Dj0bk7+KZfl7RPJ7ngaqSG+3WZVKJhFRCS8PPooxy//BM+9j5F9dBILCwqJi40hdVDXkN9usyoUzCIiEj4WLoT774err6bnQ2NZGuZLoyqie8wiIhIedu2Cq66CU06BKVPCfr3y4ajHLCIi4aF5c3j6aTj9dGja1Olq6oyCWUREQt+2bdCmDYwa5XQldU5D2SIiEpLSPV7G/PUZ3jxtIEXt2rNk2htOl1Qv1GMWEZGQk+7xMmfyXF7+z9009hdjMUz9dDM/9fZGxMzr36Ies4iIhJy0jBySNq6isb8YAwSMIfG79RGx5WZlFMwiIhJy8gsKOW37NxjAjynb2SsSttysjIayRUQk5MQd04RNLdvxZqOj2NSyHcvbJ7IqPoH4CNhyszIKZhERCS3WMvbCUxi/9wYKi3xl65Vj3C7GRsCWm5XRULaIiISOnTuhXz9S9uWSOjyR+GOPwgDxsTGkDk+M+IlfoB6ziIg4JN3jJS0jh/zgXtdjf9+FlAf+Bp9/DtaS0jO+QQTxwRTMIiJS79I9XsbPy6Kw2A+At6CQb++8Dxa/C5Mnw9lnO1yhczSULSIi9S4tI6cslAH6blnNLZ/8hw+7D4CbbnKwMucpmEVEpN4dvOxp+PpFbGrRltsH/DViD6eoKg1li4hIvYuLjcFbUEgvbzbJuVnM6H4h3/WPI/a4Fk6X5jgFs4iI1Luxg7oyZ/JcXps5DmMDFLvc/PHqiYwc1N/p0hynYBYRkXqX0jOeM7YuwB0I3mf2+7i36Q+c2gBnYR9MwSwiIvUvM5P4Tz+CqCgwhuhGjTj1qqFOVxUSFMwiIlK/duyA4cNLzld+/nlYvRr694c+fZyuLCQomEVEpH5NmAA//ABLl0Lv3nDRRU5XFFK0XEpEROrXE0/A+++XhLIcQsEsIiL1Y9ky2LsXmjeHAQOcriZkKZhFRKTurV0L558Pt9/udCUhT8EsIiJ16+ef4dJLITYWHnzQ6WpCniZ/iYhI3fH7YfRo2LoVPvkETjjB6YpCnnrMIiJSJ9I9Xt5OHgoffMD0PpeS3qS90yWFhUqD2RjTzhjzsTEm2xiz3hhza7D9AWOM1xizOvg1uNxzxhtjNhpjcowxg+ryDxARkdCT7vEyZ/JcBnk+IgBctiydOZPnku7xOl1ayKtKj9kH/N1amwAkAzcaY7oFfzfJWtsj+PU+QPB3o4BTgQuB54wxrjqoXUREQtS0Nz6n5+bVRAf8RAFuv4+em1eTlpHjdGkhr9JgttZus9auCv68G8gGfmsz02HAbGvtfmvtt8BG4MzaKFZERMLATz/x3LM30W37Jopd0fhMFMWuaJa3TzzkuEc51BFN/jLGdAR6Al8AfYGbjDHXApmU9Kp3UhLay8s9LY/fDnIREYkUPh+MGMEJe37mljOH89IZKSTnZrG8fSKr4hOIj41xusKQV+VgNsY0Bd4EbrPW7jLGTAEeBmzw+z+BPwIVnXBtK3i9McAYgPbtNSFARCQi3HknLFpE1oOTyCk+hcJiP6viEwCIcbsYO6irwwWGvirNyjbGuCkJ5RnW2nkA1trt1lq/tTYA/Jv/DVfnAe3KPb0tkH/wa1prp1prk6y1Sa1bt67J3yAiIqHglVfg6afhttvo/X+3kTo8kfjYGAwQHxtD6vBEUnSsY6Uq7TEbYwzwIpBtrX2yXHsba+224MNLgXXBn+cDM40xTwJxQBfgy1qtWkREQk+7djBiBKSlASVnLiuIj1xVhrL7AtcAWcaY1cG2e4ArjTE9KBmm3gL8GcBau94YMxfYQMmM7huttf7aLlxERJyT7vGSlpFDfkEh7Zq5uWPwqaQMHAgDBzpdWtgz1h5y+7feJSUl2czMTKfLEBGRKkj3eBk/L4vCYj+NfUXMnHUPC7r9jq6PP6Ae8hEwxqy01iYd3K6dv0RE5IikZeRQWOynV142b0wfS+/8r9jc9DitUa4l2itbRESOSH5BIb282cyZNQ53wE9xlIsdTY/VGuVaoh6ziIgckbjYGK5c/QHRgZLpQ8ZaknOziNMa5VqhYBYRkSMydlBXCpq1wGLKdvXynNhDa5RriYayRUTkiKT0jCf92Unc+HxfOn21io3dkhg5ZrgmftUSBbOIiFRNcTFccQXccAMpl1xCygu3Ol1RRNJQtoiIVM2tt8Lbb8POnU5XEtEUzCIiUrlnn4UpU+Cuu+Daa52uJqIpmEVE5Ld99FFJb/mSS+DRR52uJuIpmEVE5LdlZEBCAsyYAS6X09VEPE3+EhGRQ6R7vLw/dR4nbchkY7ckhj47lyHNmjldVoOgYBYRkQOke7y8MWkWr0y/G4OlaNls/rhvIr5mzbUkqh5oKFtERA6Q9sFX3LlgKtE2gMta3H4fPTev1l7Y9UQ9ZhEROcCFC2bTY9s3FEe5MNZS7IpmeftE7YVdTxTMIiLyP/PnM+HjF3n/5LOZdkYKyVvXsbx9IqviE4jXXtj1QsEsIiIlrIXnn+eXhNOZcMlYduJmVdtuAMS4XdoLu54omEVEpIQx8NZbHLtrF/fnFZGWkUN+QSFxsTGMHdRVE7/qiYJZRKSh+/VXuPtueOghaNECWrcmpTUKYodoVraISEPm98Po0SXbba5c6XQ1goJZRKRBSvd4GfPXZ/iyY3d4+23W3vkA/P73TpclKJhFRBqcdI+XOZPn8q+pd3Bm3np8JorUH5qS7vE6XZqgYBYRaXDSMnJI2rgKd8BX1qYNREKHgllEpIHJLyhkacfu7I9uhM9EaQOREKNZ2SIiDUleHg9+MYMHzxjFVaMeITk3SxuIhBgFs4hIQ7FrF1x8MVdu2szMbgNZFZ/AqvgEQBuIhBINZYuINATFxXDFFbB+Pe55b/KX/3ch8bExGCA+NobU4Ylatxwi1GMWEYl01sLf/gYffgjTpsEFF5CCNhAJVeoxi4hEuo0bYeZMmDABbrjB6WqkEuoxi4hEui5dYM0a6NzZ6UqkChTMIiIRJt3j5f2p8zhv2XsEYo/h6KeeJKXnSU6XJVWkYBYRiSClu3q9/J+7aewvxmK4blJvuH207imHCd1jFhGJIGkZOZyTvYzG/mIMEDCGxO/Wa1evMKJgFhGJILu3/8jgnM8AtKtXmNJQtohIBEn5fi1xu37kwQF/4ijffu3qFYYUzCIiEaTXXX/lolZd2HR0q7I27eoVXjSULSISCe6/Hz75hJSe8dz8x/O1q1cYU49ZRCTcPfkkPPQQ7NkD/fqR0jNeQRzGFMwiImEm3eMlLSOH/IJCrvnucx6a/Qhcfjk8/rjTpUktqHQo2xjTzhjzsTEm2xiz3hhza7C9hTHmI2PMN8HvxwbbjTFmsjFmozFmrTGmV13/ESIiDUW6x8v4eVkcv34Vqf+dzH1zJrKi/WnM//tEcLmcLk9qQVXuMfuAv1trE4Bk4EZjTDdgHLDQWtsFWBh8DHAR0CX4NQaYUutVi4g0UGkZOSRsWceM2RMYsfZDXDbA08kjeGzxd06XJrWk0mC21m6z1q4K/rwbyAbigWHAq8HLXgVSgj8PA16zJZYDscaYNrVeuYhIA5RfUEhybhZuv48oSjYQOf37jVqnHEGOaFa2MaYj0BP4AjjeWrsNSsIbOC54WTywtdzT8oJtB7/WGGNMpjEmc8eOHUdeuYhIA3Sqq5DzNq3AF+U6YAOROK1TjhhVnvxljGkKvAncZq3dZYw57KUVtNlDGqydCkwFSEpKOuT3IiJykF27mD7vQRr9sJkHzv8zLQp3sbx9ItkdTyNV65QjRpWC2RjjpiSUZ1hr5wWbtxtj2lhrtwWHqn8ItucB7co9vS2QX1sFi4g0SPv2QUoKsRu/Ytmkl1iyJ578gkLiYmNIHdRVy6MiSKXBbEq6xi8C2dbaJ8v9aj5wHTAx+P3tcu03GWNmA2cBv5QOeYuISOXKL4eKi41h7PknkZJ6O3z8MUyfztmjR7PU6SKlzlSlx9wXuAbIMsasDrbdQ0kgzzXG3ADkAlcEf/c+MBjYCOwF/lCrFYuIRLDS5VCFxX4AvAWFpM7+gv4bviH26adh9GiHK5S6Zqx1/vZuUlKSzczMdLoMERHH9Z24CG+5GdbGBrAmio5NXSy+90IHK5PaZoxZaa1NOrhdO3+JiISQ0mVPvbzZ3LhsLrH7dnHVqEf5bk9jhyuT+qJgFhEJIXGxMRy/fhWzZ46nUcCH30Rx2vcb+T7xkI6VRCidLiUiEkLGDurKtWv+izvgA0rWmvbdlq1jGxsQ9ZhFREJISsHX+LOXYI3Bj8EX7ab3tZfST8uhGgwFs4hIKDnmGFzn/A7uuosoj4fo/v3p16eP01VJPVIwi4iEgh9/hFatoFcvWLgQjIELNQu7IdI9ZhERp23aBImJ/ztP+fBbHksDoB6ziEg9K7+z1+lmDzNfG8vRxcUwZIjTpUkIUDCLiNSj8jt7tdj7C0/MHIfd/SOLX3qD/t26OV2ehAAFs4hIPUrLyKGw2E/vret55p3Habn3F64Z+Q+8W5to/2sBFMwiIvUqv6CQXt5sps+9j0a+YnwuF74oV9mOXyKa/CUiUo86Hu3i0nWLcPt9uLC4AgGSc7OIi41xujQJEeoxi4jUl6IiZn30T1qs/RifKxr8Popd0XhO7KGdvaSMgllEpD74fDB6NCcsWcDq8Y/y3M6mnLQhk43dkhg5Zjgp2tlLghTMIiJ1ze+HP/wB3ngDnnySHrffzlSna5KQpXvMIiJ1bdYsmD4dHnkEbr/d6WokxKnHLCJS10aPhhYtYPBgpyuRMKAes4hIXbAWHnsMNm4s2WJToSxVpGAWEalF6R4vY/76DEs79oBx48iZ+IzTJUmYUTCLiNSSdI+XOZPn8uzU2+mbuxafieIBX3vSPV6nS5MwomAWEaklaRk5/OWT6bgD/rK2nt+uJS0jx8GqJNxo8peISC3Z8dMujt/9Ez5T0ucpdkWzvH2ittuUI6JgFhGpDT4frVs2Z/g1T9Bt+2bOzFvP8vaJrIpPIF7bbcoR0FC2iEhNTZoEF1zAuHPbY49uSma7U3muzwhWxScQ43Zpu005IgpmEZGaePppuOMOaNGCS87oQOrwROJjYzBAfGwMqcMTtd2mHBENZYuIHIF0j5e0jBzyCwq5eUMGd7zzDFx6acnuXm43KT3jFcRSIwpmEZEqSvd4GT8vi4Qt67j3i3lc9M3nLOjah1/H/ZNhbrfT5UmEUDCLiFRRWkYOCVvWMWP2BBr5i/GZKF7oPZT8Rd8y7MxOTpcnEUL3mEVEqii/oJDLshbi9vtwWQvAGXnZWg4ltUrBLCJSRbeve4/Raz4gYKLwmaiydcpxWg4ltUhD2SIiVfH449zy3hQyTunLtF6XcEbeBpa3TyS742mkajmU1CIFs4hIZR55BO69F0aOZN8dE8lftJkp7U4jLjaG1EFdNQtbapWCWUTkt2Rnw/33l5yp/MorDIuOZtiZHZ2uSiKYgllE5LckJMBnn8EZZ4DL5XQ10gBo8peIyMGWLYP+/WHixJLHyckKZak36jGLiJTzyStv0/eGy4gO+PEt+YylJ3Sj3/VDnS5LGhD1mEVEgtJXbiV2wl1El56nbC0rX3uLdI/X2cKkQak0mI0xLxljfjDGrCvX9oAxxmuMWR38Glzud+ONMRuNMTnGmEF1VbiISK3y+2n8pxvonv81xVGusnXKn8WfSlpGjtPVSQNSlaHsV4B/Aa8d1D7JWvtE+QZjTDdgFHAqEAcsMMacbK3110KtIiJ1JyoKr+sonjjnapa1707y1qyy85SNdvaSelRpMFtrPzXGdKzi6w0DZltr9wPfGmM2AmcCn1e7QhGRulRYCPn50LkzLw+/Ge8v+wBY1Tah7BLt7CX1qSb3mG8yxqwNDnUfG2yLB7aWuyYv2HYIY8wYY0ymMSZzx44dNShDRKSadu+GwYOhXz/49VfGXngKMe4DZ1/HuF2M1c5eUo+qG8xTgM5AD2Ab8M9gu6ngWlvRC1hrp1prk6y1Sa1bt65mGSIi1bRzJ1xwASxZAo89BkcfTUrPeFKHJxIfG4MB4mNjSB2eqJ29pF5Va7mUtXZ76c/GmH8D7wYf5gHtyl3aFsivdnUiIrUs3eNlydOvcdcbabQo3MWqx5/nrNGjy36f0jNeQSyOqlaP2RjTptzDS4HSGdvzgVHGmMbGmE5AF+DLmpUoIlI70j1e5kyey2Ov3ctxv+4kYAxPrdut5VASUqqyXGoWJZO3uhpj8owxNwCPG2OyjDFrgfOA2wGsteuBucAG4APgRs3IFpFQkZaRQ8/Nq4GS+26uQICem1drOZSElKrMyr6yguYXf+P6R4BHalKUiEitW72aB168h5d6X0KxKxr8vrLzlPO1HEpCiLbkFJHIt2QJDBlComlCXuwJjB71CMm5/1unHK/lUBJCFMwiEtneew8uvxw6dGDNk6/x4+c/szX2BFbFl6xT1nIoCTUKZhGJGOkeL2kZOeQXFBIXG8OT0Zs4a/zfoEcPeP99BrVuTWqbA68ZO6irZmFLSDHWVrjMuF4lJSXZzMxMp8sQkTCW7vEyfl4WhcX/m2968u7tvLRpPm1f/w80a+ZgdSKHMsastNYmHdyu06VEJCKkZeRQWOynV142T777T3rlZfN1s+MZ2f8WhbKEFQ1li0hEyC8opPfW9cyZNZ5oG2BI9hJGXZWKh4TKnywSQtRjFpGI0PFoF49++CzRNgBAlA2QnJulAygk7KjHLCLhr6CA19MfoNWPuRRHuTDWUuyKxnNiD824lrCjYBaR8JeVRaucdax49F/8OzfASRsy2dgtiZFjhmvGtYQdBbOIhK9ffoFjjoFzzoEtWzijVSvOcLomkRrSPWYRCRvpHi99Jy6i07j3+PNfn2F/x07wxhslv2zVytniRGqJeswiEhZK1yknbFnHHZ73GfLVErYdczxfN47jAqeLE6lFCmYRCQtpGTkkbFnH7JnjcQd8WAwPDPgT36zfxwWXOF2dSO3RULaIhIX8gkIuXbeIRgEfBggYQ7cfNutkKIk4CmYRCQtxsTG8ddoAiqKi8ZmosiMbtU5ZIo2GskUktP38M1xzDf+4+hb+9utpjLoqtezIxuyOp5GqdcoSYdRjFpHQtXkznH02LFjAeY1/JXV4IttP7cWUPiPYfmovUocnap2yRBz1mEUkNH3xBVxyCfh8sGABnHMOKaAgloinYBaRkJHu8fL+1Hn0/TKDUWs+xNcmjgZHt6UAABoASURBVKOXZEBXDVdLw6FgFpGQkO7xMmfyXF6aPg6334c1htv7/pHBe5uS4nRxIvVI95hFJCQ8+f4Gblr0Cm6/j2gbwFjLSds2kZaR43RpIvVKwSwiziso4B/T7qZv7loCJuqA5VBapywNjYayRcRZmzbBkCGcnfsNYy+6hU0t25Uth1oVn0C81ilLA6NgFpF6ke7xkpaRQ35BIXGxMYwd1JWUwPcwaBBYy/LnZ/FubjMKi/2sik8AIMbt0nnK0uBoKFtE6lzpARTegkIs4C0oZPy8LN7b1RiSk+GLL/jdn64gdXgi8bExGCA+NkbrlKVBUo9ZROpcWkYOhcV+enmzSc5dS7P9e5n0u6t59IsfuPjdd8uuS+kZryCWBk/BLCJ1Lr+gkF7ebGbMnkBjXxFRwH6Xm6fPudrp0kRCjoayRaTOxcXGcMHXn9MkGMp+DEXRjXQAhUgF1GMWkTqXevwueqz5ACgNZTeeE3toYpdIBRTMIlLnzj3rZH7ucjL/d9oQmuXnsrFbEiPHDNf9ZJEKKJhFpFYcvBzqrgGdGJbzGVx9NXTrRou1q3jYGKfLFAl5CmYRqbHS5VCFxX4A9nvzaXvZzZC3ATp3Ljm6UaEsUiUKZhGpsfLLoYZuWMzFXy2laVEh/3flvTx09tlOlycSVhTMIlJjpcuhZs8cjzvgwwJjL7qVee2Tecjp4kTCjJZLiUiNxcXGkJybhSvgxwABE8Xxv+7UciiRalAwi0jNfP89k6K+wXNiD4qi3WUnQ2k5lEj1aChbRKrv88/h8ss5c88ern57Kbc2ieakDZlaDiVSA5UGszHmJWAI8IO19rRgWwtgDtAR2AKMsNbuNMYY4GlgMLAXuN5au6puShcRx1gLU6bAbbdB+/bw3/9y8emncXH/05yuTCTsVWUo+xXgwoPaxgELrbVdgIXBxwAXAV2CX2OAKbVTpog4Ld3jpe/ERXS6+13eOXMw3HgjXHABrFgBp5/udHkiEaPSYLbWfgr8fFDzMODV4M+vAinl2l+zJZYDscaYNrVVrIg4o3Sd8vHrV/HX5a9T4DM8c+7VpD84BY491unyRCJKde8xH2+t3QZgrd1mjDku2B4PbC13XV6wbVv1SxQRp6Vl5HD1Z3O5c8l/cAUCFLuiGT3qEbZ/9A0pvds5XZ5IRKntWdkVbe1jK7zQmDHGmExjTOaOHTtquQwRqTU+H9ekP8eExS/TyO8j2gZw+30k52aRX1DodHUiEae6wby9dIg6+P2HYHseUP6fz22B/IpewFo71VqbZK1Nat26dTXLEJE65fXCgAH85Ys3yehyFvuiG5Uth1rePlHrlEXqQHWHsucD1wETg9/fLtd+kzFmNnAW8EvpkLeIhLaDD6H4v+7NGHTdxbB3L5n/mMxt+7qQsGUdyblZLG+fSHbH00jVOmWRWmesrXCk+X8XGDML6A+0ArYD9wPpwFygPZALXGGt/Tm4XOpflMzi3gv8wVqbWVkRSUlJNjOz0stEpI4cfAgFQEx0FG/nv8fJY2+EhIRDgnvsoK5apyxSA8aYldbapEPaKwvm+qBgFnFW34mL8BYUct7GLxn/ySs8dfYo3k84l/jYGJaOG+B0eSIR6XDBrJ2/RIT8gkKuXfkuDy54HoCn3nuS75u3xkOCw5WJNDzaK1ukofP5uP/L2TwQDGUDRAUCJOdmaXKXiAMUzCIN3aRJXP/xdD456YwDZl3rEAoRZ2goW6QBqHDiVudm0Lw53HQTdOnCLx3O4Nap83QIhYjDNPlLJMIdPOO6cfF+Hlz8IkN+zKbpujXQrJnDFYo0TJr8JdJApWXkUFjsp5c3m4uzlzBw45d0/OV7Zpw7gtGNGztdnogcRMEsEuHyCwrplZfN7FnjcQd8ADw04P/x8hkpjG7UyOHqRORgmvwlEuHiYmPok7sWd8CHAfwmiia+Is24FglRCmaRCFB2VvK49+g7cRHpHm/JL957jwlntGRV556acS0SJjSULRLmDp7c5S0o5KHZK+j+yDg6vTmdwbffTtEtf+fWJtGacS0SBhTMImGudHJXqdO3fc1T7zxBh4JtcNdd8NBDpDRuTMqUmx2sUkSqSsEsEuZKz0Tu5c3mmlXvcsmGT9nerBVXjXqE2Y+Nd7g6ETlSCmaRMBcXG8Px61YxY84E3H4fGMPdF93M1tPPcro0EakGBbNIOAsEeP7HJTR/73ncfh/RNoDPRNFrx2YuH/RHp6sTkWrQrGyRcLV5M5x3Hon/fICm7dvic0XjM1H4ot30vvZSTe4SCVPqMYuEm0AAXngBxo4FlwtefpmW110Hy5fD4sVE9+9Pvz59nK5SRKpJwSwSBtI9Xt4PHjDx3cndmfjRszTr2xemTYN27Uou6tOn5EtEwpqCWSTEpXu8zHl6Dq/+526ibIDiZbO5ffg9XHj3DaS0a+t0eSJSy3SPWSTEzZj1MY++/QSNAiWTu9x+Hyd+v5m0D792ujQRqQPqMYuEKr8fJk/m1afGg7UURUWX9Jhd0Sxvn1i2fllEIouCWSRU3X47PPMMnq5nceeAv9Bm948k52axvH0iq+ITiNchFCIRScEsEkqKiuDXX+HYY+HmmyE5mR2nnEvBW+vY1rw1q+ITAIhxu3QIhUiEUjCLOKx0xvU5yz9gwHerMD16ELfov9ClC3TpQgqAMaRl5JBfUEhcbAxjB3XVOmWRCKVgFnFQusfL20+8yr9n3YfLBrBAWlRbunq8BwRvSs94BbFIA6FgFqlj6R7vYXu77057m8lv/AOXDWAAv4nCFBeRlpGjIBZpoLRcSqQOlZ6V7C0oxFJyVvL4eVmkr9wKwCpXLBtbtKPI5cZnojTjWkTUYxapSweflez2F/OHz+fSYcZq+GYVMXEncOl1T9LLm60Z1yICKJhF6lT5s5Ivy1rAOd+upv2u7fz35LNhzx7GDurK+HlZrIpP0IxrEQEUzCJ1Ki42hk6eZbzy+v1E2wAB4B/n/ZH/XjCai445hpSexwBoxrWIlFEwi9ShsYO68t0Hr+CyAQACJoqjCBzQI9aMaxEpT5O/ROrCokUwaBApXY6hx/WXs9/dWGcli0iVqMcsUkPll0P1ZBfPrJxB/IL34MQTYcsW+l0/FLp+rLOSRaRKFMwiNVC6HOrUzWtJXTaLs7auJ2Ci2PC3sXT750PQpEnJhTorWUSqSMEsUgNpGTkkbFnHf+beRxNfEQETxV9TxrK+3UCWloayiMgRUDCLVNeXXzLxhTtZd/yJuP0+DGCBk37K4yNtECIi1aRgFqnEwVtq/l/3Zgya8TTMmMGpTWP56KSzKHZFg99XtnNXnDYIEZFqUjCL/IbSe8gJW9YxLDeLDju30e++T/BHGVz33MOyi67l9YxvWXdC57Kdu7I7nkaqNggRkWpSMIv8hrSMHBK+XceMORNw+30ALGt/Ok+NvIt5j1zJEMB3dFPSMhoxJT6BuNgYUrVBiIjUQI2C2RizBdgN+AGftTbJGNMCmAN0BLYAI6y1O2tWpogDrOXUFR/zyAfP0MhXjAuLz0SxvMPpeEzzssu0QYiI1Kba2GDkPGttD2ttUvDxOGChtbYLsDD4WCS8fPYZ/O53TJ33D/a73Phc0Qec/qR7yCJSV+piKHsY0D/486vAYuDuOngfkboxejTMnAlt2uCZMJFrA6fSZetXuocsIvXCWGur/2RjvgV2UrJK5AVr7VRjTIG1NrbcNTuttcdW8NwxwBiA9u3b9/7uu++qXYdIdaV7vLw/dR49PZ/i6Xkug8cMJ2XxXNi/H265BY466pBZ2TpkQkRqgzFmZbnR5jI17TH3tdbmG2OOAz4yxnxV1Sdaa6cCUwGSkpKq/68DkWpK93jJmDiN515/CJcNULwinev2+eCWEQcEr+4hi0h9qtE9ZmttfvD7D8BbwJnAdmNMG4Dg9x9qWqRIrdu+ncJbbmPyGw/jsgEMEGUD9Ny8mrSMHKerE5EGrNrBbIw52hjTrPRn4AJgHTAfuC542XXA2zUtUqRWWQsDBzJi6Tw+7dizZHJXuYld+dq1S0QcVJOh7OOBt4wxpa8z01r7gTFmBTDXGHMDkAtcUfMyRaqn9B7yaWuW0qF4N2by0wztcxI8+yyj529hubsVvbzZZRO7VsUnEK8Z1yLioGoHs7V2M9C9gvafgIE1KUqkNqR7vLz3+Es8N/dBogN+AB6/L45A2gOk9OvHqOYnsWZeFqviE1gVnwBAjNvFWM24FhEH1cY6ZpHQU1RE4S238ezcB3EH/BjAb6IwxUVl95BTesaTOjyR+NgYDBAfG0Pq8ERN9BIRR2lLToksv/wCxxwDbjedt2SztEN3+uRmER3wV3gPWTOuRSTUKJglrJXeQz4rcyHdduWT5M3GvXkTHHccd/5lErm7i3QPWUTCioJZwla6x0vmfWlMef8ZomwAgPe7nYtZt53BA47jjsHdGK97yCISZhTMErIq23HrtdmfMu+9p7GAAXwmivWtOvL2lzsYPICya7Vrl4iEEwWzhKTSc5ALi0tmU3sLCnlwzgo6vv4aPfb/CP/8Jx7TnNR+13Pb0pm4/T7dQxaRiKBglpCUlpFDYbGfXt5sBmz8kta/7uT3G7+kReEuSE6GoiLiYmN4IflyVrQ7VfeQRSRiKJglJOUXFNLLm83sWeNx+30ArIjvxhP9rmXu9LvAGMYO6qp7yCIScRTMElqKiuDNN7n45220z80iKhAou3+8uHMS3sQzoGS3Od1DFpGIpGAWx5QudTppQyY7OpzMDTE/ccr8WfD999x90aXcdeI5FC+bDcH7x54TexzSG9Y9ZBGJNApmcUS6x8ucyXN5afo4GvuKMJTMrP7+dwM44aWXaDdoECPXbOPWJtGctCGTjd2SGDlmuEJYRCKeglnqzGGXO/34I9/d9yhnFBTg9vuIAgLAv8+4lFcvvomlFw0Agr3hKTc7+jeIiNQ37ZUtdaJ0uZO3oBALeHfu5a0np5N3YQrEx3Pre8+xL7oxxa5ofCaK/dGN+KDr2TpyUUQaPPWYpU6UX+503qYVXLZuIXG7f2JXk6bwl79wtTmdz5qcQGa7blrqJCJSjoJZqu2wQ9X79tFj+Uf80buBq1Z/ULbc6amzR/F88hV89c/hXO7xslJLnUREDqFglmo5ZGeunXuZ8/Qcuu/x0GnhuzxbUMBPMc1x+31E2wA+E0VRdCNato4FtNRJRORwFMxSLeWHqpNzs2hb8D1Xrf2QQncTGHUFS88ezLOf5/HizAmHXe6kpU4iIodSMMuR83oZ9NFsRq7NoOPObbgCfnxRLib3GcnUsy5j3aQr6AvsOMvLrUc11nInEZEjoGCWCpXf/GNjtySGXH8xQ1d+ALNnw2ef8X/Wsv3oY4n2+3BhIQD73I055viWZa+h5U4iIkdOwSyHKN384+Xp43D7iilaNpsb9u7ngvcep8lxreCBB1iQ2J8X56/kpenjyk52qmhnLhEROTIK5gbqsDOqN28m974nePyzeWU7crn9PnpsyeLyPz/Hu/+4DIzhfGBPx87amUtEpJYpmBug0hnVCVvWMSw3i+XtEhn/axEJzz5O1xcncwuwJfYEfFEujLVl5xyv98eUHSABGqoWEakLCuYGaPI7axi99E3GLX4Flw2w3+XmqisfZdLRnXn+ySe5YltrVkQdWzbjWpt/iIjUHwVzhDp48tbgMcNJiS2Cv/yF/y78mMb+YiwlB0dEB/wk52Yxpc8IuP1iRnu8rNPmHyIijlAwR6B0j5f5T7zKlDn34/b78H02k+v2+XD9v4u55PvveSt5KF81bsHdn75aNnFreftE4oI9Ym3+ISLiHAVzpHn4YbpMncm0vByisABEB3z03LyaiUu6c8maNTTxeJkzL4u1cSeXDVVndzyNVG3+ISLiOAVzGEr3ePnouTn87ssPadooisS2x9DxrVklv1yxgr02itcTB5Ky4RNcAX9Zj7j05Kb/9YgbMSU+gbjYGFLVIxYRCQnGWut0DSQlJdnMzEynywgtn38OixdD//7Qp09Zc+Yjz3DMU09w0o+5GMACG1u1Z8N7ixl2Ziewlr6PfYy3oLDCyVtLxw1w6A8SEZHyjDErrbVJB7erx1wHDrtGuIo+ffFNzvrLlUT7fGAMRS1bE5P5BXTowIJVWxhdtK9s4pbfRPFWt/68vejbkmA2hrGDujJek7dERMKSgrmWHXLqUkEh4+dlAf8bQi4/Y/rbU3oy4rxunNf/dDjhBJZOmcXZN15NtA0AYK0lzzThu5VbOL9DB17oMoAvj2rDjNkTDpi4VTpMXf59NHlLRCT8KJhrWempS+UVFvtJy8ghpWc8736Ww8933sNzX6bjsgH49DXMVFh710Oc/th9pOVGcf0p5zA4ZylRNkCxK5q7B93E9q/9nA/ExcawigRGj3rkN9cYa/KWiEh4anDBXNNh5sqen19QSLP9v3JhzlLO25RJsSuaYwt3s6xjdxg3gKcWbmTBF/PKhqIDwLunnMO0qJN5G1hjm3Lb0LG85h1yQPCaYI9Yw9QiIpEtooK5stCs6jDz4V6j/FaWl3/rYesxx7N4xUJObuOnW2In+POfiYuN4b/3j6B50V6gZHLW5mPjiTn6KAA2FcINw+/jX/MfKxuKfjlpKGttU6CkR+wtKDwgeEvby9epYWoRkcgUMcFcldAtHWY+eLZy6TBzWfB+m8WVm1bw3bFtWJTZjM4doklsewxp+7uTsGUdr8+4C9fBs9kHDYI//5mxg7qy9LXeXJD9GS4sfhPF2z1+z4kP3wOUBOzCLmcddii6tEdcfjj84B6xhqlFRCJXxATzb4Zut1ZQUEB+cAnRzFn34Pb7CERF8W7X3xFwuWD1VNJ6jCkJ3pnjSu7/lnfSSeRf9hTDcrMI7tuBH8N/el1MWr/rWP/k5UBJaH5y150Uj/kS6yvGF+2m97WX0i8YpJUNRatHLCLSsEVMMJeG7ozZE2jkK8YABU2acpRvP4wvAqD9wx+QnJtFI38xUYAr4OeSr5bwY/NWcHQn8jvsDQZvSfL6MbzW+2Ke+t3VrJk0grjHPmZ5+0SKot1lw9Dzu/Uj9rgWB9TS7/qh0PVjWLyY6P796VduHXJVglc9YhGRhitigjkuNobkz7Nw+324sASA3NgTWH9yT0YP6g6xsfy9dxdmbepB0bJZRPtLdsT64+hURt46kpSe8cRNXHRI8L6T0I+mJ7T+3/rgX4sOGIY+eCvLMn36HLAxSHkKXhEROZw6C2ZjzIXA04ALmGatnVhX7wUlQ8RzcnpQvGw2BEP1sQv/yshbRkAwBIcCgUYjuKVJdNmpSyPHDC8LycqCV1tZiohIXauTLTmNMS7ga+D3QB6wArjSWruhoutra0vOCo86PMLQrOlyKhERkao43JacdRXMfYAHrLWDgo/HA1hrUyu6Xntli4hIQ3O4YI6qo/eLB7aWe5wXbCtf0BhjTKYxJnPHjh11VIaIiEh4qatgNhW0HdA1t9ZOtdYmWWuTWrduXUdliIiIhJe6CuY8oF25x22B/Dp6LxERkYhRV8G8AuhijOlkjGkEjALm19F7iYiIRIw6WS5lrfUZY24CMihZLvWStXZ9XbyXiIhIJKmzdczW2veB9+vq9UVERCJRXQ1li4iISDUomEVEREKIgllERCSEKJhFRERCiIJZREQkhNTJXtlHXIQxO4DvavElWwE/1uLrNVT6HGtOn2HN6TOsOX2GNVcXn2EHa+0hW1+GRDDXNmNMZkUbg8uR0edYc/oMa06fYc3pM6y5+vwMNZQtIiISQhTMIiIiISRSg3mq0wVECH2ONafPsOb0GdacPsOaq7fPMCLvMYuIiISrSO0xi4iIhKWIC2ZjzIXGmBxjzEZjzDin6wk3xph2xpiPjTHZxpj1xphbna4pXBljXMYYjzHmXadrCVfGmFhjzBvGmK+C/5vs43RN4cYYc3vwv+V1xphZxpgmTtcU6owxLxljfjDGrCvX1sIY85Ex5pvg92Pr6v0jKpiNMS7gWeAioBtwpTGmm7NVhR0f8HdrbQKQDNyoz7DabgWynS4izD0NfGCtPQXojj7PI2KMiQduAZKstadRcgzvKGerCguvABce1DYOWGit7QIsDD6uExEVzMCZwEZr7WZrbREwGxjmcE1hxVq7zVq7Kvjzbkr+jzDe2arCjzGmLXAxMM3pWsKVMaY5cC7wIoC1tshaW+BsVWEpGogxxkQDRwH5DtcT8qy1nwI/H9Q8DHg1+POrQEpdvX+kBXM8sLXc4zwUKtVmjOkI9AS+cLaSsPQUcBcQcLqQMHYisAN4OXhLYJox5miniwon1lov8ASQC2wDfrHWfuhsVWHreGvtNijpwADH1dUbRVowmwraNO28GowxTYE3gdustbucriecGGOGAD9Ya1c6XUuYiwZ6AVOstT2BX6nD4cNIFLwPOgzoBMQBRxtjrna2KqlMpAVzHtCu3OO2aNjmiBlj3JSE8gxr7Tyn6wlDfYGhxpgtlNxOGWCMme5sSWEpD8iz1paO2LxBSVBL1Z0PfGut3WGtLQbmAWc7XFO42m6MaQMQ/P5DXb1RpAXzCqCLMaaTMaYRJZMc5jtcU1gxxhhK7ullW2ufdLqecGStHW+tbWut7UjJ/wYXWWvVSzlC1trvga3GmK7BpoHABgdLCke5QLIx5qjgf9sD0QS66poPXBf8+Trg7bp6o+i6emEnWGt9xpibgAxKZh++ZK1d73BZ4aYvcA2QZYxZHWy7x1r7voM1ScN1MzAj+A/tzcAfHK4nrFhrvzDGvAGsomTFhQftAlYpY8wsoD/QyhiTB9wPTATmGmNuoOQfPFfU2ftr5y8REZHQEWlD2SIiImFNwSwiIhJCFMwiIiIhRMEsIiISQhTMIiIiIUTBLCIiEkIUzCIiIiFEwSwiIhJC/j+XF34oVq7VcwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_fitted = results.fittedvalues\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(x,y,'o',label='data')  # 原始数据\n",
    "ax.plot(x,y_fitted,'r--.',label='OLS')  #拟合数据\n",
    "ax.legend(loc='best')\n",
    "plt.show()  # 可以看到图已经非常精确了"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**分类变量：**\n",
    "假设分类变量有3个取值（a,b,c），比如考试成绩有3个等级。a就是（1,0,0），b（0,1,0），c（0,0,1），这个时候就需要3个系数β0，β1，β2，也就是β0x0+β1x1+β2x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
       "       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
       "       0, 0, 0, 0, 0, 0])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nsample = 50\n",
    "groups = np.zeros(nsample, int)\n",
    "groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [1., 0., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 0., 1.]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#亚变量，转化成类onehot的形式，变成0，1的形式\n",
    "groups[20:40] = 1\n",
    "groups[40:] = 2\n",
    "dummy = sm.categorical(groups, drop=True)\n",
    "dummy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.994</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.994</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   2589.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 10 Nov 2020</td> <th>  Prob (F-statistic):</th> <td>2.80e-51</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>21:53:52</td>     <th>  Log-Likelihood:    </th> <td> -74.545</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    50</td>      <th>  AIC:               </th> <td>   157.1</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    46</td>      <th>  BIC:               </th> <td>   164.7</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    7.9672</td> <td>    0.635</td> <td>   12.551</td> <td> 0.000</td> <td>    6.689</td> <td>    9.245</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    2.0457</td> <td>    0.073</td> <td>   28.010</td> <td> 0.000</td> <td>    1.899</td> <td>    2.193</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>    <td>   -0.0252</td> <td>    0.403</td> <td>   -0.063</td> <td> 0.950</td> <td>   -0.835</td> <td>    0.785</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x3</th>    <td>    2.4251</td> <td>    0.336</td> <td>    7.209</td> <td> 0.000</td> <td>    1.748</td> <td>    3.102</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x4</th>    <td>    5.5673</td> <td>    0.758</td> <td>    7.346</td> <td> 0.000</td> <td>    4.042</td> <td>    7.093</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 0.042</td> <th>  Durbin-Watson:     </th> <td>   1.779</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.979</td> <th>  Jarque-Bera (JB):  </th> <td>   0.224</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.028</td> <th>  Prob(JB):          </th> <td>   0.894</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.677</td> <th>  Cond. No.          </th> <td>1.28e+17</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The smallest eigenvalue is 4.16e-31. This might indicate that there are<br/>strong multicollinearity problems or that the design matrix is singular."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.994\n",
       "Model:                            OLS   Adj. R-squared:                  0.994\n",
       "Method:                 Least Squares   F-statistic:                     2589.\n",
       "Date:                Tue, 10 Nov 2020   Prob (F-statistic):           2.80e-51\n",
       "Time:                        21:53:52   Log-Likelihood:                -74.545\n",
       "No. Observations:                  50   AIC:                             157.1\n",
       "Df Residuals:                      46   BIC:                             164.7\n",
       "Df Model:                           3                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          7.9672      0.635     12.551      0.000       6.689       9.245\n",
       "x1             2.0457      0.073     28.010      0.000       1.899       2.193\n",
       "x2            -0.0252      0.403     -0.063      0.950      -0.835       0.785\n",
       "x3             2.4251      0.336      7.209      0.000       1.748       3.102\n",
       "x4             5.5673      0.758      7.346      0.000       4.042       7.093\n",
       "==============================================================================\n",
       "Omnibus:                        0.042   Durbin-Watson:                   1.779\n",
       "Prob(Omnibus):                  0.979   Jarque-Bera (JB):                0.224\n",
       "Skew:                          -0.028   Prob(JB):                        0.894\n",
       "Kurtosis:                       2.677   Cond. No.                     1.28e+17\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The smallest eigenvalue is 4.16e-31. This might indicate that there are\n",
       "strong multicollinearity problems or that the design matrix is singular.\n",
       "\"\"\""
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Y = 5+2X+2Z1+6*Z2+9*Z3，Z1、Z2、Z3分别表示取a、b、c\n",
    "x = np.linspace(0,20,nsample)\n",
    "X = np.column_stack((x, dummy))\n",
    "X = sm.add_constant(X)\n",
    "beta = [5,2,3,6,9]  #即上面的假设\n",
    "e = np.random.normal(size=nsample)\n",
    "y = np.dot(X, beta) + e\n",
    "result = sm.OLS(y,X).fit()\n",
    "result.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "        <script type=\"text/javascript\">\n",
       "        window.PlotlyConfig = {MathJaxConfig: 'local'};\n",
       "        if (window.MathJax) {MathJax.Hub.Config({SVG: {font: \"STIX-Web\"}});}\n",
       "        if (typeof require !== 'undefined') {\n",
       "        require.undef(\"plotly\");\n",
       "        requirejs.config({\n",
       "            paths: {\n",
       "                'plotly': ['https://cdn.plot.ly/plotly-latest.min']\n",
       "            }\n",
       "        });\n",
       "        require(['plotly'], function(Plotly) {\n",
       "            window._Plotly = Plotly;\n",
       "        });\n",
       "        }\n",
       "        </script>\n",
       "        "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from plotly.offline import init_notebook_mode,iplot  # 可以交互的包，鼠标对着有结果\n",
    "import plotly.graph_objs as go\n",
    "\n",
    "init_notebook_mode(connected=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
